#pragma once

#include "vector.h"

namespace g2 {

template <typename T>
class basic_line_segment;

typedef basic_line_segment<double> line_segment;
typedef basic_line_segment<float> line_segmentf;
typedef basic_line_segment<long double> line_segmentl;

template <typename T>
class basic_line_segment {
public:
	typedef basic_line_segment<T> line_segment;
	typedef basic_vector<T> vector;
	vector v0, v1;
public:
	line_segment()
		{ }
	line_segment(vector const & start, vector const & finish)
		: v0(start), v1(finish) { }
	template <class U> 
	line_segment(basic_line_segment<U> const & b) 
		: v0(vector(b.v0)), v1(vector(b.v1)) { }
public:
	vector center() const
		{ return (v0 + v1) / T(2); }
	vector get_vector() const
		{ return v1 - v0; }
	T square_length() const
		{ return sqr(v1 - v0); }
	T length() const
		{ return absolute(v1 - v0); }
	/// returns [squared] distance from a point to the segment
	T square_distance(vector const & p) const;
	T distance(vector const & p) const
		{ return sqrt(square_distance(p)); }
	/// returns [squared] distance from a segment to the segment
	T square_distance(line_segment const & b) const;
	T distance(line_segment const & b) const
		{ return sqrt(square_distance(b)); }
};

template <typename T>
inline T basic_line_segment<T>::square_distance(vector const & p) const
{
	T s = scalar_product(v1 - v0, p - v0);
	if (s < 0)
		return sqr(p - v0);
	T sqlen = square_length();
	if (s > sqlen)
		return sqr(p - v1);
	return sqr(outer_product(v1 - v0, p - v0)) / sqlen;
}

template <typename T>
inline T basic_line_segment<T>::square_distance(line_segment const & b) const
{
	// Consider substitution:
	// e0 := b.v1 - b.v0,	
	// e1 :=   v1 -   v0.
	// The intersection point of lines, containing segments, satisfies equations:
	// [e0, p - b.v0] = 0,
	// [e1, p -   v0] = 0.
	// We can rewrite the equations using scalar products:
	// (orth(e0), p - b.v0) = 0,
	// (orth(e1), p -   v0) = 0.
	// Since biorthogonal basis for (orth(e0), orth(e1)) is (e1/[e0, e1], -e0/[e0, e1]),
	// any vector x can be decomposed
	// x = -(orth(e1), x)/[e0, e1]*e0 + (orth(e0), x)/[e0, e1]*e1.
	// Find the value of (e0, p - b.v0), bearing in mind that (orth(e0), p - b.v0) = 0,
	// (e0, p - b.v0) = (e0, p - v0 + v0 - b.v0) = 
	//  = -([e1, p - v0] + [e1, v0 - b.v0]) / [e0, e1]*(e0)^2 = 
	//  = [e1, b.v0 - v0]) /[e0, e1]*(e0)^2.
	// The segments do intersect each other if and only if they are not parallel ([e0, e1] != 0) and
	// 0 <= (e0, p - b.v0) <= (e0)^2, therefore
	// practical conditions of intersection are
	// [e0, e1] != 0 && [e1, b.v0 - v0]*[e0, e1] >= 0 && [e1, b.v0 - v0]^2 <= [e0, e1]^2

	vector e0 = b.v1 - b.v0;
	vector e1 =   v1 -   v0;
	T x = outer_product(e1, b.v0 - v0);
	T ee = outer_product(e0, e1);
	if (ee != 0 && x * ee >= 0 && sqr(x) <= sqr(ee))
		return 0;

	return 
		std::min(square_distance(b.v0),
			std::min(square_distance(b.v1),
				std::min(b.square_distance(v0),	b.square_distance(v1))));
}

} // g2
